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Theories in which gravity is weaker on cosmological scales have been proposed 

^O . to explain the observed acceleration of the universe. The nonlinear regime in such 

O; 

^^ I theories is not well studied, though it is likely that observational tests of structure 

^ ' formation will lie in this regime. A class of alternative gravity theories may be 

o: 

^^ ■ approximated by modifying Poisson's equation. We have run N-body simulations of 

''^ ! a set of such models to study the nonlinear clustering of matter on 1-100 Mpc scales. 

_^ I We find that nonlinear gravity enhances the deviations of the power spectrum of these 

> ■ 

QQ ' models from standard gravity. This occurs due to mode-coupling, so that models 

(^ , with an excess or deficit of large-scale power (at k < 0.2 Mpc^^) lead to deviations 

^; 

^^ in the power spectrum at smaller scales as well (up to /c ~ 1 Mpc ), even though 

^^ ' the linear spectra match very closely on the smaller scales. This makes it easier to 

o ■ distinguish such models from general relativity using the three-dimensional power 

2 ! spectrum probed by galaxy surveys and the weak lensing power spectrum. If the 

52 I potential for light deflection is modifled in the same way as the potential that affects 

^ ■ the dark matter, then weak lensing constrains deviations from gravity even more 

strongly. 

Our simulations show that even with a modifled potential, gravitational evolution 
is approximately universal, and our conclusions extend to models with modiflcations 
on scales of 1-20 Mpc. Based on this, the Peacock-Dodds approach can be adapted 
to get an analytical fit for the nonlinear power spectra of alternative gravity models, 
though the recent Smith et al formula is less successful. We also use a way of 
measuring projected power spectra from simulations that lowers the sample variance, 
so that fewer realizations are needed to reach a desired level of accuracy. 
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I. INTRODUCTION 



Observations of Type la supernovae, along with observations of the cosmic microwave 



background and_ 
is accelerating 



arge-scale structure, have established that the expansion of the universe 

n 

|3|]. Einstein's theory of gravity, and a cosmological model that in- 
cludes dark matter, baryons and radiation, cannot explain this cosmic acceleration. The 
explanation may involve the existence of an exotic form of energy density, or the breakdown 
of general relativity (GR) on large scales. These explanations for cosmic acceleration are 
known as "dark energy" and "alternative gravity" approaches, respectively. 

The rate at which the universe expands is predicted by the Friedmann equation H^ = 
87tGp/3 (for a spatially fiat universe), which is derived from the Einstein equations and 
the metric. In order to reproduce the cosmic acceleration, we have to modify the Einstein 
equation: 

Rf„y - -Rgf,^ = SttGT^^. (1) 

The left hand side describes the curvature of spacetime due to gravity while the right hand 
side describes its sources. Given an observed expansion history that one wishes to describe, 
an alternative gravity (AG) theory will attempt to explain it via a modification of the left- 
hand side while a dark energy (DE) theory introduces a new term on the right-hand side that 
gives the desired acceleration. While one could imagine two different kinds of modifications 
that gave the same expansion history, they will have different effects on the growth of 
structure by gravitational clustering: a smooth dark energy (DE) affects the growth of 
structure only by changing the expansion rate, while AG affects it by direct modification 
of the gravitational interaction. The linear regime growth factor G{a) is scale independent 
in DE models; an AG modification will produce a different solution, which in general is 



scale-dependent: G{k, a) as described in Sec. |TT] (though see Jacobs et al. 4] for weak scale 
dependence in standard gravity due to the effect of inhomogeneities). 

There are several potential ways of modifying GR: adding nonlinear terms in the Ricci 
scalar R to the gravitational action, coupling i? to a scalar field as in Brans-Dicke gravity, 
having gravity operate in a higher dimensional universe on large scales as suggested by brane 
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cosmology, introducing scalar and vector degrees of freedom and so on. The tensor-vector- 
scalar theory of Bekenstein [5| , the ghost condensate theory of Arkani-Hamed et al. |6| , and 
the five-dimensional theory of Dvali, Gabadadze, and Porrati [7, DGP gravity] have drawn 
recent attention. Our interest is in the class of theories such as DGP that aim to reproduce 
cosmic acceleration by weakening gravity on large-scales. The theory is then likely to have 
testable consequences on scales probed by large-scale structure. 

Lue et al. [8| have derived the linear growth of perturbations in DGP gravity. Further, 
Lue et al. |9(] argued that generic gravity theories that obey Birkhoff' s theorem and mimic 
cosmic acceleration lead to the suppression of the growth of large-scale density perturbations 
at the level of ~ 5%, similar to DGP. Predictions for the exact linear growth as a function 
of scale, or for scale-dependent growth in the nonlinear regime, do not yet exist in such 
theories. 



Newton's Law has been directly measured from millimeter to solar system scales [10|, lll |. 
To constrain possible deviations from Eq. ([T]) on cosmological scales requires geometric infor- 
mation (distance measurements to objects of known redshift) and information on the growth 
of structure. Since the observed cosmic acceleration occurs at low redshift, observational 
constraints at z < 1 are needed to learn about its origin. Geometric information at low red- 
shift has been measured most cleanly by the measurement of the luminosity distances dilz) 
to type la supernovae. Weak lensing (WL) can measure low-redshift distance information as 



well as the evolution of the growth factor, especially with tomography 12|. Baryon acous- 
tic oscillations in the power spectrum of galaxies or other probes can measure the angular 
diameter distance at low redshifts. The CMB (cosmic microwave background) measures 
geometric information (the angular diameter distance to last scattering at z = 1089), the 



shape and amplitude of the primordial power spectrum, and matter /radiation content 13|. 
Thus the CMB anchors the cosmological model at high redshift, and by comparing to it. 
Type la SN, weak lensing, baryon oscillations and galaxy cluster measurements constrain the 
effects of dark energy or AG. Currently, combining the CMB, SN^nd large-scale structure 
information has led to a best-fit "standard" cosmological model [2|. This model is usually 
described in terms of standard gravity (GR) and dark energy; current constraints on the 
equation of state of dark energy are consistent with a cosmological constant, but a possible 
time evolving dark energy is not well constrained. 

Some studies have considered the constraints on alternative gravity from existing and 



planned survey data [e.g. 



14, 


15. 


16 



17 



mm 



m l2l|, 122, l23|, 124 . Lue et al. (8| and some of 



these authors have examined the hnear regime growth factor in DGP gravity and computed 
its consequences for low-redshift galaxy power spectra and weak lensing observables. 

In this paper we study the consequences of modifying Poisson's equation for the growth of 
structure in the nonlinear regime. The model we consider for the modified Poisson equation 
may be regarded as an approximate description of a complete AG theory over a range of 
scales (though not all AG theories will be described by our approach to gravitational clus- 
tering). While our AG model is not derivable from a covariant, consistent theory of gravity, 
it has the merit that we can use N-body simulations to study the the full nonlinear evolution 
of structure. We are interested in the amplitude and scale dependence of modified growth 
over the range 1-100 Mpc, where cosmological observations can probe gravity effectively. An 
AG theory that provides an expansion history similar to the ACDM accelerating universe 
is likely to alter Newton's Law on these scales. We note that a similar modification of the 
gravitational potential was considered by White and Kochanek 25l]. They followed a some- 
what different approach to lensing (changing the deflection angle relation while retaining 
the growth of structure as in standard gravity) and calculated the consequences on smaller 
scales (0.01-10 Mpc), as their focus was on modifications that produce fiat rotation curves 
for galaxies without the need for dark matter. 

In the linear regime, analytic calculations can explore the effects of a modification of 
Poisson's equation on the growth of structure. Recent efforts in this direction were made by 
Shirata et al. 2m\ and Sealfon et al. 27|]. However observations have significant information 
in the small scale nonlinear regime, so it is necessary to develop predictions for this regime. 
Moreover, we know from perturbation theory that quasilinear effects propagate power from 
large to small scales, so altering gravity on large scales is likely to affect smaller structure as 
well. To obtain accurate nonlinear predictions, we use N-body simulations to determine the 
effect of modifications to Poisson's equation in the non-linear regime on structure formation. 
We quantify the extent to which such a modification would be constrained by galaxy and 
WL surveys. 

In Sec. ITT] we describe the formalism that describes the growth of perturbations due to 
gravity. Sec. IIIII contains details on our numerical simulations and predictions for three- 
dimensional and lensing power spectra. We describe our results in Sec. HV] and conclude in 
SeclVl 



II. LINEAR REGIME 

In Eulerian coordinates the equations that govern the behavior of mass fluctuations are 
given by recasting the fluid equations in expanding coordinates, or simply by conservation 
of stress-energy V^T^'' = 0. If we linearize these equations, the resulting second order 
differential equation describes the growth of the fractional overdensity 6{r, t), or equivalently, 
its Fourier transform 6{k,t): 

5 + 2H6 = ^, (2) 

~ ~ k"^ ~ 

6 + 2H5 = — ^0, (3) 

where a{t) is the expansion scale factor and gives the Hubble parameter as H{t) = a/a. The 
Fourier transformed Poisson equation in comoving coordinates reads 

where k is the comoving wavenumber and is the gravitational potential. In this work our 
continuous Fourier transform conventions are 

^(k) = I dh5{v)e'^\ (5) 

i0J(k)e---, (6) 

A DE modiflcation will change Eq. ([3]) via the time derivatives and the Hubble parameter 
H = a/a on the left hand side; for such a model one separates Eq. ([3]) by letting (5(k, t) = 
6(k)G{t) and then solving for the growth factor G{t). An AG modiflcation will change the 
potential on the right-hand side via Eq. (jl]). We can see from this that if an AG modiflcation 
is made, i.e. if fc^0(k, t) ~ (5(k, t)/(k, t) for some non-trivial /(k, t), Eq. ([3]) will no longer 
be separable, and hence the growth factor will become scale dependent, i.e. one must allow 
^(k,t) = 5(k)G(k,t). Eq. (^ then becomes 



3 Hq iZmO 



G(k, t) + 2HG{k, t) = -^L-^G(k, t)/(k, t). (7) 



If /(k, t) -^ 1 then ^(k, t) —>■ G{t) as in standard gravity. We note that a purely time- 
dependent modiflcation f{t) (such as a time- dependent Newton's constant), can be accom- 
modated without a scale-dependent growth factor. 
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FIG. 1: Ratios of linear theory alternative gravity (AG) power spectra to standard gravity at z = 0, 
for different parameterizations of Eq. Q. The dotted lines are the matter power spectrum ratios 
-P5,ait/-P<5,std = {Gaitik,t) /Gstd{t))'^ ■ The dashed lines are the ratios of the velocity power spectra 
Pv,ait/Pv,std = iGa.itik,t)/Gstdit))'^- The expansion history is fixed to be the same as for a ACDM 
universe. 

A. Modified Poisson Equation 



A theory of gravity that makes gravrty weaker or stronger over a^arrge of lerrgth scales can 

be approximated by modifying Poisson's equation. Sealfon et al. [27] investigated the effect 
of a Yukawa-type (adding an exponential term) and power-law modifications to the Poisson 
equation, solving for the scale-dependent growth factor G{\i.,t) under the assumption that 
the modification was a small perturbation to the Newtonian potential. Shirata et al. 26 ] 
followed by obtaining G(k, t) for a Yukawa modification of arbitrary strength. They adopted 
the Peacock-Dodds (PD) prescription [23] — an approach we check with our simulations — 
to obtain the non-linear matter power spectrum from the linear solution. With a prescription 
for galaxy biasing, this enabled them to predict galaxy power spectra, which they used with 
their AG model to constrain model SDSS galaxy power spectrum measurements 29|. 



Both Shirata et al. 26|] and Sealfon et al. 27|] consider the real-space potential 



!>ait r 



G / d'^Y 



3 . P(r') 



1 + a 1 -e 



In Fourier space, this becomes 



v^0a.(k) = y^ 



<5(k) 



1 + a 



l + (|k|r,/a)2j 



(8) 



(9) 



Note that this will result in a scale dependent growth factor when plugged into Eq. ([2]) 
above. In Fig. [H we can see the effect of this modification on the linear theory power 
spectrum; there we have plotted the ratio of the AG matter and velocity linear power 
spectra to the corresponding standard linear spectrum at redshift 2; = for a few different 
parameterizations of Eq. ([9]). Throughout this work, we fix the background expansion to be 
the same as ACDM, i.e. we take 



H\a) = H^{n,noa-'' + n 



AOj 



(10) 



In this study we do not allow the acceleration to vary with a or r^; we let it be fixed solely by 
A, because an AG theory would need to predict an accelerating expansion not too different 
from that given by ACDM in order to fit the supernova data. For Fig. [Hand our simulations 
we take Qj^o = 0.3, and Q^o = 0.7. 

The matter power spectrum Ps{k) oc G'^{k, t), while the velocity power spectrum Pv{k) oc 
G^{k,t), so the ratios of the growth factors give us the ratios of the linear spectra starting 
from the same initial spectrum. We solve for the growth factor using Eq. (j7]). The velocity 
power spectra show a more pronounced difference than the matter spectra (a factor of 2-5 
larger deviation at A; ~ 0.05 Mpc^^). It may be worth exploring the measurement of large- 
scale peculiar velocities via the kinetic SZ effect or distance measurements on galaxies to 
test gravity. 

On large and small scales, we get the limiting behavior 



r > rs, 0alt ^ (1 + a)0newton, 
r <^ Ts, 0alt ^ 0newton- 



(11) 

(12) 



So for positive a, the alternative gravity pot ential is strong er than Newtonian gravity on 
large scales and unchanged on small scales. 



Shirata et al. 



take Tg to be a fixed physical 



length, that is, in comoving units it changes with redshift; consequently, at early times 



8 

when a <^ 1, Tg becomes very large. At our simulation starting point oi z = 50, Vg is much 
larger than the boxsize, and so the linear spectra are virtually identical. Hence both the 
alternative and the standard gravity simulations start from the same initial conditions. We 
examine the effect of the initial conditions further in Sec. IIVI 

We note that this modification cannot extend to arbitrarily large scales or early times. 
We regard it as an approximate description of gravity on length scales well below the horizon 
at low redshifts. 

III. NUMERICAL SIMULATIONS 

While Eq. ([2]) is useful for describing linear-regime density fluctuations on large scales, if 
we want a more complete description we have to turn to numerical simulations. An N-body 
code simulates the evolution of structure by evolving a large number of particles interacting 
by gravity. These are evolved from early times (2; ^ 1) up to the present, with particle 
positions and velocities outputted at regular intervals. 

An efficient N-Body solver must compute the forces on a large number of particles simul- 
taneously so that the equations of motion can be integrated forward in time. We use a basic 
particle-mesh (PM) solver for this purpose, which interpolates the particles onto a grid and 
then computes the potential via Fourier transform. More advanced techniques (e.g. P^M 
and tree codes) are available, which provide larger dynamic range in exchange for greater 
complexity and computation time. We use PM simulations to simulate modified gravity in 
the quasilinear to moderately nonlinear regime where observations can test models without 
needing to consider astrophysical/baryonic effects. Due to the lower CPU costs of PM sim- 
ulations, we were able to run a large number of realizations to reduce sample variance on 
the power spectra, which would have been prohibitive with the other methods. Our code 
is based on the PM code of Klypin and Holtzman [30|], which was designed and tested for 



DM simulations 



31 



32j, and kindly made available by A. Klypin. We set up the initial 
conditions by displacing particles from a regular grid using a realization of the linear power 
spectrum. 



A. Discrete Poisson Equation 

In the PM simulation, the equations of motion are discretized on a grid, starting with 
the Poisson equation. Since we modify this equation for the AG simulations, we give the 
explicit formulae here. We define the second derivative operator in 1-D as 

We define the unnormalized discrete Fourier transform as 

N-l 



r=0 
N-l 



i2-Krk/N 



fc=0 

Combining the previous two equations leads to the discrete Fourier space expression 

2Tik 



[V20], = 0fe X 2 



cos ■ 



So the 3-D discrete Poisson equation for standard gravity in Fourier space reads 

7 _ 3 //q ^mO ^k 



(13) 



where 



Gi 



(14) 



cos — — + cos — — + cos — — - 3 . (15) 
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a Gk 




2nk^ 


2nky 2nk^ 
N +^°^ N 
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B. Particle-Mesh Simulation Parameters 

The parameters which determine the dynamic range of a PM simulation are the boxsize 
Lbox, the number of particles N^, and the Fourier grid size Ag. The grid spacing should 
be smaller than the mean particle spacing to preserve the small scale resolution; a common 
choice which we adopt is to take Ag = 8Ap, i.e. twice the number of particles per dimension. 
Lbox must be large enough so that we have enough power in the linear regime to get accurate 
power spectra, however increasing Lbox decreases one's ability to resolve small-scale structure 
for fixed Ap. Our computational resources fixed Ag = 256'^ and Ap = 128^; we chose 
-^box = lOO/i"^ ^ 140 Mpc for our boxsize. The wavenumber corresponding to Lbox is 
/Cjiiin ~ 0.04 Mpc~^; the comoving distance to the sources is x(z = 1) ^ 3.3 Gpc, so the 
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FIG. 2: Diinensionless 3D power (A^(A;) = k^Ps{k)/2iT'^ , left panel) and convergence power 
(A^(^) = £'^P^{£)/2Tr, right panel) from N-body simulation with a standard (Newtonian) grav- 
ity model. Also shown are the Smith et al. (solid black line) and Peacock-Dodds (solid red line) 
fitting formulae. The points with error bars are an average over 8 realizations. Three sets of 
simulations are shown with varying resolution due to differences in the total number of particles; 
the total number of grid points in each set is A'^g = 8Np. 
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FIG. 3: Dimensionless 3D power and convergence power spectra from simulation, shown together 
with the Smith and PD fitting formulae as in Fig. [2l Three different boxsizes are shown: lOO/i^^, 
200/i-\ and 500/1"^ Mpc. 
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angular wavenumber is i^i^ = k^mX ~ 145, corresponding to a field of view ~ 2.5° on a 
side. 

To check our results in the non-linear regime, and to insure that we chose Lbox large 
enough, we tested our prediction for the 3D power for standard gravity against the Smith 



et al. 33| fitting formula. We find in Fig. [2] that for our runs of A^p = 128^ particles on a 
A^g = 256^ grid, our results are limited by resolution at physical scales of about 1.0 Mpc, and 
angular scales of £ ~ 1000. The power spectrum A^(/c) in Fig. ([2]) is identical to the linear 
theory prediction on large scales; we have about a decade of power in the linear regime at 
z = 0. With our sources at Zg = I, the weak lensing weight function W{x) iii Eq. ( TT9l) peaks 
at z ^ 0.4 where the linear regime extends to smaller scales. So we can be confident that 
for purposes of measuring the lensing power spectrum our choice of Lbox = I00h~^ Mpc is 
large enough. 

Fig. [2] and Fig. [3] show how the resolution limit behaves as we change the number of 
particles or the boxsize in the simulation. We note that since we are not including direct 
particle-particle effects (as in P^M-type codes) we don't need to be concerned with explicit 
force softening, and that the shot-noise contribution is very small on the scales that we 
resolve. 

The particle mass in our simulation rrip is given by 

mp = fi^0Pcr,0 (^] (16) 

For our simulations with A^p = 128^ and Lbox = lOOh^^ Mpc, we get nip = 1.1 x 10^° Mq. 

We ran all of our simulations in a ACDM background cosmology. They were started at 
redshift z = 50. We used as = 1.0, Hq = 70, ^a = 0.7, and Qm = 0.3 as our cosmological 
parameters. We ran 8 realizations for our runs with standard gravity and for each of the 
alternative gravity models. 

C. Convergence Po^ver Spectrum 



The output of N-Body simulations 



maps of cosmological weak lensing 3J, 



las been used to make model shear and convergence 



351]. Unlike the angular power spectrum of the CMB, 
which can be computed analytically, weak lensing involves the deflection of light by large 
scale structure at low redshift, which has already undergone significant nonlinear collapse, 
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FIG. 4: Comparison of the dimensionless convergence power using our Limber approximation 
method (red triangles) and the standard multiple lens plane method (blue squares), each averaged 
over 8 realizations. Computing P^ in the standard method (using Eq. (j22|) ) uses only the 2D 
modes of the density field at each redshift slice. We included the full 3D modes by using the 
Limber equation for P^ (Eq. (j24p l. The scatter is significantly smaller, as shown by the error bars 



28| 3D power 



on the red symbols. The black curve is computed using the Peacock and Dodds 
spectrum. 

so model WL maps and power spectra must be computed numerically. 

In the standard multiple lens plane algorithm, the convergence and shear maps are com- 
puted by filling the light cone from the observer to the source galaxies with matter from 
N-body simulation outputs. Starting at the source redshift, one can output three orthogonal 
2D-projections (slices) of the density field at length intervals equal to the box size. By pick- 
ing a projection at each redshift and randomly translating the boxes, we tile approximately 
uncorrelated mass distributions along a light cone for each realization. The box size at the 
source redshift determines the field of view simulated. 



Following Bartelmann and Schneider 36|], in the thin-lens approximation for a flat uni- 
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verse, the convergence k{6) is a sum over slices of the transverse gradient of the potential 
at comoving distances Xi (with sources fixed at Xs)'- 

0) = \Y1 Ll,o.W{x^, Xs)Vl<P{x^O, X.)- (17) 



K 

C 

I 



In standard general relativity, we can use the Poisson equation, the real-space version of 
Eq. (jlj), to obtain 

Kstd{0) = — ^ — 2^ Lbo^W{xi, Xs) , . , (18) 

where Lbox is the tile boxsize, and 

W{x.,Xs)^ ^'^^'~^'^ (19) 

Xs 

is the weak lensing weight function. If the photons feel the same modified potential as the 
dark matter, then we must use a modified Poisson equation such as Eq. ([H]) instead; note 
that this is not true in some theories, e.g. in the fifth force-type modification of Nusser 



et al. [37|]. The converse has also been considered: White and Kochanek [25| constrained 
the parameters of a model where the potential affecting photons was modified, but structure 
was assumed to have formed as in GR. Since we are interested in scales {£ > 100), we use 
the flat-sky approximation and define the 2D Fourier transform as 

k{£) = [d^0K{9)e'^-'', (20) 



dH . 



<^) = / TTTT^^We" (21) 



(27r)2 

Assuming the mass distribution in different redshift slices is uncorrelated, the convergence 
power spectrum in this approximation is a sum over the 2D power spectra of the projected 
densities: 

Here we have included the function /(k, a), introduced in Eq. ([7]), that describes the mod- 
ification of gravity in the Poisson equation; for the model we consider in this paper (i.e. 
Eq. O), 

We can already see that if an AG modification changes the potential felt by dark matter 
and photons in the same way, then WL results can provide a stronger constraint on AG. 
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In Eq. ( !22|) . the modification will affect the two-point function (|k(^)P) twice: once via the 
modified growth of structure from the {\Sf^{i)\'^) term, and once via the /^(fc, a) term that 
comes directly from the modification of the potential. We show in Sec. HV] that for our 
model, these effects of AG combine to produce a dramatic effect. 

Computing the convergence power spectrum from Eq. (l22l) uses only the 2D modes of 
the density field at each redshift slice; a large part of the information in each simulation box 
is lost by the projection. In order to reduce the scatter in the simulations, we included the 
full 3D modes in our computation of P^i^) from our data. We accomplished this by using 
the Limber approximation to compute P^ directly as an integral over the 3D matter power 
spectrum Ps: 

P.,„ . £^AxE ^/^ (* - ^."(X.)) Ps U - ^.X.) . (24) 

4c^ i Xia^iXi) \ Xi J \ Xi J 

Traditionally, Eq. ( l24l) is used to compute -Pk(^) when one has an estimate of Ps{k) from 
the halo model or other fitting formula. Instead we use Eq. ( I24l) with Ps measured from the 
PM simulations at redshifts Zi. This results in a significant gain in our signal-to- noise for 
Pk as shown in Fig. HI the error bars in the standard method based on measuring P^ from 
K, maps are significantly larger than in our method. 

For each realization we obtain the power spectrum in 100h~^ Mpc boxes tiled along the 
line of sight. Using the discrete version of Equation (12^ we get Pk(£) for each realization 
from Ps{k) binned in spherical shells in wavenumber k; the bins we need at each redshift 
are given hj k = £/xi- Eight realizations give us our estimate of the scatter in our results. 
We take our results to be valid up to wavenumber k ^ 1 Mpc~^, which is about k^yq/2; for 
Pk this corresponds to £ ~ 1000. Additionally, we have checked the validity of our modified 
gravity simulations for larger boxsizes. Since at the start of our simulation, the comoving 
scale rs{z = 50) ~ 250 Mpc is larger than Lbox = 100h~^ Mpc, but at 2 = 0, r^ ^ i^box, we 
needed to insure that no numerical artifacts are produced when r^ crosses the box scale. In 
Fig. [5], we find that this is indeed the case. 

IV. SIMULATION RESULTS 

We have run ensembles of simulations using standard gravity and the AG potential given 
by Eqs. (jH [9]) with five different sets of parameters: a = 0.5, rg = 5 Mpc; a = — 0.5,rs = 
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FIG. 5: Dimensionless 3D power and convergence power spectra from simulation. The colored 
pentagons are simulation data produced using a 400h~^ Mpc box, the squares have a 200h^^ Mpc 
box, and the triangles have lOOh^^ Mpc (which we use in the main portion of our simulations). 
The simulations are consistent with each other up to their resolution limit, regardless of whether 
or not rg is initially outside the box. 



5 Mpc; a = 0.5, rg = 10 Mpc; and a = 
are within the 2a range of constraints set by 



0.5, Ts = 10 M pc. These sets of parameters 



Shir at a et al. 



using SDSS data; the last two 



models given have the smallest deviation from standard gravity, for these models the linear 
spectrum differs by 20% at A; = 0.05 Mpc^^. We also consider a model that has significantly 
less power on large scales: a = —1.0, r^ = 5 Mpc. 

The left panel of Fig. E] shows the 3D power for standard gravity and the five AG models, 
while the left panel of Fig. [7] shows them as ratios to standard gravity. There is a statistically 
significant difference between the models at the smallest scales resolved by our simulations, 
where the general trend is that models with larger |a| and smaller r^ are more different 
from standard gravity. The overall shape of the linear and nonlinear AG spectra remains 
similar to the shape of the standard gravity spectra; as expected, the positive a models 
have excess power on large scales, while the models with negative a have less large-scale 
power compared to standard gravity. The nonlinear scale of each of the models is around 
k = 0.1-0.2 Mpc~^, similar to standard gravity; on scales smaller than this, the linear 
spectra asymptote to the standard gravity linear curve. However, the 3D nonlinear spectra 
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FIG. 6: Dimensionless 3D power (A^(/c) = k^P{k)/27r'^) and convergence power (A^(£) = 
l'^Pn{t)/2'K) for standard and alternative gravity models. The black curves for standard grav- 
ity are computed using the Peacock and Dodds [28] fitting formula for the nonlinear 3D power, 
and using it in the Limber integral for P^. The dotted curves are predictions from linear theory. 
The solid curves are our analytic fits for the nonlinear spectra of modified gravity models (as indi- 
cated by the legends), while the symbols show simulation measurements. The two subpanels for A^ 
illustrate two choices for how WL power spectra are affected by a modified potential: the curves 
in the left subpanel are computed assuming a GR deflection law, i.e. setting /^ = 1 in Eq. (j24p . 
whereas those on the right are computed with a modified potential for photons. The error bars 
on the solid curves in the right panels show the expected statistical errors from a future lensing 
survey (see text in Section 4 for details). 



show clearly that changing gravity on large scales propagates into the nonlinear regime: in 
the nonlinear region k = 0.5-1.0 Mpc^^, where the linear spectra are within a few percent 
of standard gravity, the nonlinear 3D spectra differ by 10% or more. 

The Pk plots in Figs. [^|71 are split up into two subpanels: in the left subpanel, photons 
are not affected by the modification of the potential (but the growth of structure is still 
altered), corresponding to setting /^ = 1 in Eq. fl2^ . while in the right panel photons feel 
the modified potential. Each P^ is given by the limber integral over the corresponding Ps 
in the left panel, the difference is only whether the /^ term in Eq. flMl) is included. If light 
deflection and structure formation are both affected by the modified potential, then the 
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FIG. 7: Ratio plots of dimensionless 3D power (A^(/c) = k^P{k)/2'K'^) and convergence power 
(A^(£) = i'^Pfi(£)/2ir). The solid curves show the ratio of the Peacock-Dodds prediction for AG with 
standard gravity. The dotted curves are the linear theory ratios. We divided our AG simulation 
power spectra by the corresponding points from the standard gravity simulation to obtain the 
points. As in Fig. [6l the left P^ subpanel shows the WL result for the GR deflection law, while 
the right subpanel shows that the deviation of AG from the standard gravity result is substantially 
enhanced if the deflection law is modified along with the growth of structure. The error bars on the 
horizontal line in the right panels show the expected statistical errors from future lensing surveys 
(see text for details). 



difference between the P^ subpanels in Figs. [^|71 shows that an AG model would be much 
more easily ruled out. Because the modified potential for photons and the modified growth 
of structure can reinforce each other, WL statistics potentially have more power to constrain 
AG theories than large-scale structure observations (such as the galaxy power spectrum) do: 
in the left panel of Fig. [71 we see that the cyan and green models {\a\ = 0.5, r^ = 10 Mpc) 
are perhaps just barely detectable using the power spectrum; but if light deflection as well 
as large-scale structure is modified, these models are easily ruled out. 

For Pfj as well the differences from standard gravity extend to smaller angular scales; at 
our resolution limit of £ ~ 1000, there are still observable differences at the level of several 
standard deviations for our fiducial survey. The error bars on the black lines in each of the 
right panels are the statistical errors from a hypothetical lensing survey with /sky = 0.1 and 
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FIG. 8: Comparing the effect of changing the initial conditions to changing the potential. The ICs 
at 2; = 50 are on the left, and the z = outputs are on the right. The cyan and magenta points are 
simulations that have a modified initial spectrum shape but are evolved with the standard gravity 
potential, while the red and blue points are evolved using the AG potential with the correct initial 
spectrum (as in the rest of the paper). The two pairs of points for each model (blue and cyan 
points, and the red and magenta points), are within errors of each other at the end of the simulation 
(the z = right panel). Hence there is an approximate degeneracy between the shape of the initial 
power spectrum and the shape of the potential during its evolution under gravity. 

^gai = 40 arcmin^^. The general trends are the same as for Pg: the models with larger |a| 
and smaller r^ exhibit larger deviations from standard gravity, and the linear and nonlinear 
spectra are more different on large angular scales. As is evident from Fig. [6l while the linear 
Pk asymptotes to the the standard gravity result on small scales, in the observable regime 
{i > 100), Pk has nonlinear contributions. Here we have not performed a full parameter 
analysis that would involve CMB priors and variations of all relevant parameters. 

A. Analytical Approximations 



From the power spectra in Fig. [H] and Fig. [71 one cannot tell whether the observed 
differences on small scales at late times are a result of the changed non-linear evolution in 
the AG models, or whether the structure formed under the influence of normal gravity and 
merely started from different initial conditions. 
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FIG. 9: Dimensionless power for negative and positive a models at different redshift, compared 

n n 

to the Smith et al. 331] (blue, dashed) and Peacock and Dodds [2^ (red, solid) fitting formulae. 
The z = outputs have been translated by a factor of 4 in the y-direction for legibility. For 
negative a, PD fits better than Smith, while for positive a, the simulation points lie between the 
two predictions. 

To answer this question we ran simulations whose initial conditions had the same shape 
as the late-time AG linear power spectrum, but were evolved using the standard gravity 
potential. The results, shown in Fig. [HI are striking: the 3D power spectra of simulations, 
which had the z = linear shape of the AG model at the initial time, came out the same as 
the regular AG runs. Recall that since r^ is a fixed physical length scale in the AG model, 
at the start of the simulation aX z = 50, it is 250 — 500Mpc in comoving coordinates. This 
is larger than the simulation box size of 100/i~^ Mpc. So we would expect the AG power 
spectrum at z = 50 to be identical to that of standard gravity. 

The plots in Fig. [8] show a kind of universality in cold dark matter structure formation: 
the way non-linear structures form is not uniquely determined by specifying the detailed 
shape of the potential. Our simulations show that the DM power spectrum by itself cannot 
distinguish between changing the shape of the initial power spectrum and changing the 
shape of the gravitational potential. 

This is further revealed by testing the power spectra measured in simulations against an- 
aljd;ical fitting functions that have been calibrated for standard gravity. The Peacock-Dodds 
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(PD) formula 28|, using a mapping of length scales between the linear and nonlinear regime, 



gives the dimensionless nonlinear 3D power spectrum A^ as a function of the dimensionless 
inear power A^l to around 10% accuracy (when compared to simulations). Shirata et al. 



26( 1 use the PD formula to extend their results to non-linear scales. 

We have tested the PD formula with AG simulations by replacing the standard gravity 
linear power spectrum in the formula with the linear spectrum from the AG model. The 
results for Ps are shown in the left panel of Fig. [6|, in which the dashed lines are produced 
using the Peacock-Dodds formula. The right panel shows P^ fits generated by integrating 
the 3D power along the line of sight using the Limber approximation. The fits in Fig. [6] are 
accurate at about the 10% level, with the accuracy apparently better for models closer to 
standard gravity. From the ratio plots in Fig. [TJ we also see that the positive a model (more 
large scale power) does worse. 



The fitting formula of lSmith et al.l performs the same task as the PD formula in that, given 
a wavenumber k and redshift z, it provides an estimate of the non-linear power spectrum 
/\1{k,z). The Smith formula is inspired by the halo model: it breaks up the non-linear 
power into two pieces, a quasi-linear term and a one-halo like term. We adapted the Smith 
formula for use with alternative gravity by using the AG linear spectrum, but we see in 
Fig. [9] that PD is a much better fit to the data for negative a (for the positive-a model that 
we have tested, the Smith and PD fits are comparable at low redshift). This may be due to 
the use of two separate terms in the Smith formula, one of which is calibrated on the basis 
of the shapes of CDM halos in standard gravity. It was designed and tested for simulations 
with scale-free initial spectra or CDM initial conditions, but not for initial spectra with 
a very different shape, such as those produced by a scale-dependent growth factor, initial 
conditions with a shape like those in Fig. [HI 

To summarize, our results suggest that the nonlinear power spectrum in alternative grav- 
ity models is captured completely by the change in the linear growth factor. This result is 
consistent with the approach used in the PD fitting function. A similar result is shown in the 



recent study of hinder and White [38[ , who found that nonlinear spectra for a class of dark 
energy models can be accurately described by appropriate choice of length and time scales. 
We must emphasize, however, that our results are only valid for a certain range of scales. 
We can estimate over what range our results should be valid by examining what happens to 
the comoving scale r^ and the nonlinear scale as the simulation progresses. The comoving 
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wavenumber corresponding to r^ is defined by ks{z) = 2tt/{{1 + z)rs). For a r^ = 5 Mpc 
model, Fig. [TOl shows how the wavenumber kg ranges from approximately 1.3 Mpc^^ at the 
present to 0.16 Mpc^^ at z = 6.78. Inspection also reveals that the nonlinear wavenumber 
fcni (the scale where the linear and nonlinear spectra begin to diverge) ranges from approx- 
imately fcni = 0.2 Mpc~^ at the present to fcni = 0.5 Mpc~^ at z = 6.78 for the a = —1, 
Tg = 5 Mpc model. So we see that r^ actually starts out larger than r^i, but ends up well 
inside the nonlinear scale as time goes by; the redshift where they coincide depends on the 
values of r^ and alpha chosen, for a = — 1, r^ = 5 Mpc one can see that it occurs at about 
z ^ 3.2, whereas for the a = 0.5, r^ = 10 Mpc model, the scales cross later at a redshift 
of around z ~ 1.9. As long as the boxsize and resolution of our simulations can adequately 
capture the dynamics of r^ and r^i over a range of redshifts, i.e. for any model where the 
nonlinear scale and r^ cross well before redshift z = 0, the conclusions we reach from our 
simulation should be valid; we estimate that this would be true for a range of r^ at least 
from about 1 Mpc to 20 Mpc. 



It is interesting that quasilinear perturbation theory [39|, |40|] would have suggested some 
departures from this universality: equations for the second order density and velocity fields 
contain the linear fields as well as the V0 term in the Euler equation. So the dependence 
of the second order terms on the scale dependent function f{k,t) introduced in Eq. ([7]) 
would not be completely determined by the linear solution. Our results show that this 
quasilinear departure is likely very small. Generally speaking, the halo model description, 
assuming it describes alternative gravity models, is consistent with universality: halo bias, 
the linear spectrum and halo mass function depend only on the linear growth and initial 
power spectrum. Since the small scales where halo structure may play a role are not well 
probed by our simulations, it is not surprising that we find an approximate universality in 



the non 
ACDM 



inear power spectrum. The PD formula gives results close to the halo model for 



33|. 



V. DISCUSSION 

We have performed N-body simulations of large scale structure formation with a modified 
Newtonian potential in a ACDM background. This is intended to approximate alternative 
gravity theories that are designed to match the observed acceleration of the universe. We 
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FIG. 10: Dimensionless 3D power spectrum {A'j{k, z) = k^P{k, z)/27r'^) as a function of wavenum- 
ber and redshift for two different modified gravity models. The solid lines are the linear spectra 
while the points and error bars are the average of eight realizations from simulation; the vertical 
long-dashed lines indicate the value of ks{z) = 2TTa{z)/rs, which is the comoving wavenumber 
corresponding to the physical scale r^. The nonlinear scale k^\{z) can be read off the graph by 
observing where the linear and nonlinear power spectra begin to diverge. We can see that in each 
model, ks < k^i initially, but at z = 0, kg > k^i in each case. 



focus on the quasilinear and non-linear regime of clustering at low redshift. Our simulations 
resolve the 3D power spectrum of matter on scales oi k ^ 0.05-1.0 Mpc~^. We used the 
3D simulations to compute the weak lensing power spectrum over angular wavenumbers 
i ^ 100-1000. We used a technique for this that reduces the scatter in measurements from 
simulations (described in Section 3.3). The range of scales we studied is expected to be 
observable with high accuracy with planned surveys. The nonlinear modification of the 
power spectra ranges from 10% effects in the quasilinear regime to an order of magnitude at 
the small scale end. While the accuracy of our simulated spectra is typically 10% over this 
range, the relative accuracy for predictions of different gravity models is significantly better. 
We find that nonlinear effects propagate the difference between the power spectra to scales 
where the linear spectra are nearly identical. This is the expected effect of mode coupling 
in nonlinear gravitational evolution 391]. The result is that at scales of fc ~ 0.5 Mpc^^, 
the modified gravity power spectra differ by over 10%, while the linear spectra are within 
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5%. Similar differences are seen in the weak lensing spectra at £ ~ 500, and if the potential 
for photons is modified in a similar way as the potential for matter, the effects reinforce 
each other and the differences in the weak lensing spectra are much larger. These scales are 
of great interest because the observational errors are expected to be small and theoretical 
interpretation can be made without modeling of non-gravitational effects (at least for the 
lensing spectra). We compare the differences between models to the expected statistical 
errors from a wide area lensing survey to show that it should be possible to directly constrain 
the parameters of an alternative gravity scenario. Models within the 2-cr limits of current 
galaxy surveys would be easily distinguished by future lensing measurements. A detailed 
study of this is left for future work. 

Our results are consistent with a universality in nonlinear gravity, which makes the non- 
linear power spectrum a function only of the initial (Gaussian) conditions and linear growth, 
for modifications on length scales for which there are likely to be precise observations in the 
near future. We find that the lensing and 3D power spectra cannot distinguish between sim- 
ulations which were started with appropriately modified initial conditions and evolved with 
standard gravity, and those which were evolved with a modified gravitational potential. This 
means that for observationally accessible scales, there is a degeneracy between the shape 
of the potential used to evolve the simulation and the shape of the initial conditions. For 
simple modified gravity models, it means that another constraint on the primordial power 
spectrum (such as the CMB), or measurements at multiple redshifts, must be used to test 
for modified gravity. It may also be that the structure of dark matter halos, not well probed 
by our simulations, are different for AG models. This is an interesting topic for future work 
with higher resolution simulations. The skewness or bispectrum may also help distinguish 
alternative gravity models, as suggested by Bernardeau [4l|] and Sealfon et al. 27|, since its 
dependence on the scale dependent function f{k,t) in Eq. (jTj) is different. 

We tested analytical approximations to the nonlinear spectra for modified gravity models. 
We found that while the Peacock-Dodds fitting formula was accurate to 10-20% in compar- 
ison to the simulated spectra, its relative accuracy between different models is significantly 
better and within the errors of our measurements. The Smith et al. formula does somewhat 
worse for the models studied. 

Our simulations are a useful first step in studying the effects that an alternative gravity 
model has on large scale structure formation. It would be of great interest to simulate 
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the dynamics of a full alternative gravity theory, instead of our approach (which may be 
a convenient approximation but is not even covariant). However, a full alternative gravity 
theory that one can simulate is hard to come by; e.g. going beyond the linear regime for the 



DGP model 



A 



even in the quasilinear regime, is an unsolved problem. 
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